import keras
from keras import models
import shap
import lightgbm as lgb
import sklearn
import numpy as np

# model_path = "tongji-1.0-1-5.h5"
# model = keras.models.load_model(model_path)
# explainer = shap.GradientExplainer(model)

X,y = shap.datasets.adult()
X_display,y_display = shap.datasets.adult(display=True)

X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(X, y, test_size=0.2, random_state=7)

model = lgb.LGBMClassifier(random_state=1234)
model.fit(X_train, y_train)

explainer = shap.TreeExplainer(model)
# 最新版本的shap对于LGBMClassifier得到的shap_values为两个数组的列表，即两个分类的结果，这里使用分类1的结果
shap_values = explainer.shap_values(X)
if len(shap_values) == 2:
    shap_values = shap_values[1]
# 如果报错，更新lightgbm
shap_interaction_values = explainer.shap_interaction_values(X)
print(shap_interaction_values)
print('解释模型的常数:', explainer.expected_value)
print('训练样本预测值的log odds ratio的均值:', np.log(model.predict_proba(X_train)[:, 1]/ (1 - model.predict_proba(X_train)[:, 1])).mean())
print('常数与归因值之和:', explainer.expected_value + shap_values[0].sum())
print('预测值:', model.predict_proba(X.iloc[0:1])[:, 1][0])
print('预测值的log odds ratio:', np.log(model.predict_proba(X.iloc[0:1])[:, 1][0] / (1 - model.predict_proba(X.iloc[0:1])[:, 1][0])))

# 输出如下：
# 解释模型的常数: -2.359295492472404
# 训练样本预测值的log odds ratio的均值: -2.3592954924724046
# 常数与归因值之和: -4.342723316250527
# 预测值: 0.012834215330737286
# 预测值的log odds ratio: -4.342723316250522

print('第一个特征的shap values:', shap_values[0][0])
print('第一个特征的shap interaction values:', shap_interaction_values[0][0])
print('第一个特征的shap interaction values和:', shap_interaction_values[0][0].sum())

# 输出如下：
# 第一个特征的shap values: 0.6792772188209449
# 第一个特征的shap interaction values: [ 0.70680724 -0.05856777  0.07723929  0.00210557 -0.04619057  0.07169773  0.01072035  0.0353209  -0.09792147 -0.00315917 -0.01261797 -0.00615691]
# 第一个特征的shap interaction values和: 0.6792772188209456

'''
force plot是针对单个样本预测的解释，它可以将shap values可视化为force，
每个特征值都是一个增加或减少预测的force，预测从基线开始，基线是解释模型的常数，
每个归因值是一个箭头，增加（正值）或减少（负值）预测。
'''
shap.initjs()
shap.force_plot(explainer.expected_value, shap_values[0], X.iloc[0,:])

'''
summary plot是针对全部样本预测的解释，
有两种图，一种是取每个特征的shap values的平均绝对值来获得标准条形图，这个其实就是全局重要度，
另一种是通过散点简单绘制每个样本的每个特征的shap values，通过颜色可以看到特征值大小与预测影响之间的关系，同时展示其特征值分布。
'''
shap.summary_plot(shap_values, X, plot_type="bar")
shap.summary_plot(shap_values, X)

'''
如果要看特征值大小与预测影响之间的关系使用dependence plot更合适，
dependence plot清楚地展示了单个特征是如何影响模型的预测结果的，
dependence plot同样有多种使用方式，一种是查看某个特征是如何影响到模型预测结果的，
另一种是一个特征是如何和另一个特征交互影响到模型预测结果的。
'''
shap.dependence_plot('Age', shap_values, X, interaction_index=None)
# 分别为散点图的x与y即可绘制得到，运行下面代码，可以得到与上面一摸一样的图：
# plt.figure(figsize=(7.5, 5))
# plt.scatter(X['Age'], shap_values[:, 0], s=10, alpha=1)

#接下来看一个特征是如何和另一个特征交互影响到模型预测结果的，这里以Age和Education-Num为例子，运行下面代码：
shap.dependence_plot('Age', shap_values, X, interaction_index='Education-Num')

# 使散点图的垂直没那么分散，使用下面代码可以绘制：
shap.dependence_plot(('Age', 'Age'), shap_interaction_values, X, interaction_index=None)
# 这个图的绘制也很简单， 因为shap interaction values的维度是N×M×M，
# 我们只需用[:, feature_index, feature_index]切片就可以得到​，用下面代码可以得到上图的结果：
# plt.figure(figsize=(7.5, 5))
# plt.scatter(X['Age'], shap_interaction_values[:, 0, 0], s=10, alpha=1)

'''
并不是所有模型都支持类别型变量的，需要将特征进行one-hot编码，这时候就需要自己去处理shap_values和shap_interaction_values，
除非将类别下变量进行label编码当作数值型去处理，然后再绘制归因关系时特征值使用原始值。
'''